Differentially expressed genes in orbital adipose/connective tissue of thyroid-associated orbitopathy

Background Thyroid-associated orbitopathy (TAO) is a disease associated with autoimmune thyroid disorders and it can lead to proptosis, diplopia, and vision-threatening compressive optic neuropathy. To comprehensively understand the molecular mechanisms underlying orbital adipogenesis in TAO, we characterize the intrinsic molecular properties of orbital adipose/connective tissue from patients with TAO and control individuals. Methods RNA sequencing analysis (RNA-seq) was performed to measure the gene expression of orbital adipose/connective tissues of TAO patients. Differentially expressed genes (DEGs) were detected and analyzed through Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, and Gene Set Enrichment Analysis (GSEA). The protein–protein interaction (PPI) network was constructed using the STRING database, and hub genes were identified by the Cytoscape plug-in, cytoHubba. We validated several top DEGs through quantitative real-time polymerase chain reaction (qRT–PCR). Results We identified 183 DEGs in adipose tissue between TAO patients (n = 3) and control patients (n = 3) through RNA sequencing, including 114 upregulated genes and 69 downregulated genes. The PPI network of these DEGs had 202 nodes and 743 edges. PCR-based validation results of orbital adipose tissue showed multiple top-ranked genes in TAO patients (n = 4) are immune and inflammatory response genes compared with the control individual (n = 4). They include ceruloplasmin isoform x3 (CP), alkaline tissue-nonspecific isozyme isoform x1 (ALPL), and angiotensinogen (AGT), which were overrepresented by 2.27- to 6.40-fold. Meanwhile, protein mab-21-like 1 (MAB21L1), phosphoinositide 3-kinase gamma-subunit (PIK3C2G), and clavesin-2 (CLVS2) decreased by 2.6% to 32.8%. R-spondin 1 (RSPO1), which is related to oogonia differentiation and developmental angiogenesis, was significantly downregulated in the orbital muscle tissues of patients with TAO compared with the control groups (P = 0.024). Conclusions Our results suggest that there are genetic differences in orbital adipose-connective tissues derived from TAO patients. The upregulation of the inflammatory response in orbital fat of TAO may be consistent with the clinical phenotype like eyelid edema, exophthalmos, and excess tearing. Downregulation of MAB21L1, PIK3C2G, and CLVS2 in TAO tissue demonstrates dysregulation of differentiation, oxidative stress, and developmental pathways.


INTRODUCTION
Thyroid-associated orbitopathy (TAO), also called Graves' ophthalmopathy, is a category of autoimmune diseases associated with thyroid dysfunction (Bahn, 2010).A prominent feature of TAO is the expansion of orbital tissue, comprising both extraocular adipose and muscle tissues (Garrity & Bahn, 2006).The swollen soft tissues are the result of the accumulation of nonsulfated glycosaminoglycan, inflammation, hyaluronan, and the activation of local fibroblasts (Berchner-Pfannschmidt et al., 2016).If left untreated, the expansion of orbital tissue can result in orbital congestion, significant exophthalmos, compressive neuropathy, and even lead to vision loss causing a serious decline in quality of life (Wang et al., 2021).In the last several decades, rehabilitative orbital decompression surgery has been the standard treatment for the stable stage of TAO.This surgical approach aims to mitigate proptosis, alleviate orbital congestion, and enhance the aesthetic appearance of the orbital region.Consequently, it serves as a means to ameliorate the quality of life for individuals afflicted with TAO (Bartalena, 2013).
The activation of orbital fibroblasts plays a key role in the immune process of TAO pathogenesis (Naik et al., 2010).Under pathological conditions, orbital fibroblasts will express functional molecules, such as thyrotropin receptor, the receptor of insulin-like growth factor, and CD40, and continue to differentiate into adipocytes and myofibroblasts closely related to disease progression.Most of the current studies focus on isolating and establishing primary orbital fibroblasts and conducting further immune research related to various pathological mechanisms of TAO (Hammond et al., 2021;Jang et al., 2019).However, limited research has been conducted concerning the direct detection of gene expression within the orbital adipose/connective tissue of TAO patients utilizing highthroughput sequencing methods.This issue emphasizes the importance of comprehending the underlying mechanism(s) of orbital adipogenesis to identify therapeutic approaches for the prevention or treatment of TAO.
The transcriptome refers to the sum of all RNA transcripts for a specific tissue or cell in a certain developmental state or functional condition, including messenger RNA (mRNA), noncoding RNAs, and small RNAs.Screening the specific genes that play a key role in disease among many differentially expressed genes (DEGs) has become a key research goal (Chen et al., 2021).Bioinformatics analysis based on gene expression profiles may screen hub genes and regulatory pathways, which play an important role in the early diagnosis of TAO and the establishment of early warning mechanisms (Kim et al., 2021).
In this study, DEGs were identified based on high-throughput RNA sequencing data of tissues from TAO and control subjects to explore the pathogenesis of TAO.Then, Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Set Enrichment Analysis (GSEA) pathway analyses were obtained to predict the functions of these DEGs.The expression patterns of some DEGs were confirmed by qRT-PCR.

Subjects and tissue samples
All human studies were conducted according to the Declaration of Helsinki principles and were approved by the Ethics Committee of the Affiliated Wuxi People's Hospital of Nanjing Medical University (identifier, KY23013).We collected human orbital adipose/connective tissues from 43 to 80-year-old patients with TAO undergoing routine resection of prolapsed orbital fat in the Department of Ophthalmology, the Affiliated Wuxi People's Hospital of Nanjing Medical University, from July 2021 to August 2022.The demographics of the patients are presented in Table 1 and Fig. S1.All TAO patients included in this study were diagnosed according to Bartley's criteria, and tissues of control individuals obtained in plastic surgery were collected as control samples.All patients provided written informed consent.

Bulk RNA sequencing analysis (RNA-Seq)
The total RNA in tissues were extracted.To ensure the quality of the samples for transcriptome sequencing, the concentration and integrity of RNA samples were checked using a Nanodrop ND-2000 spectrophotometer and an Agilent Bioanalyzer 2100/4200, respectively.The qualified RNA samples were used for mRNA preparation and cDNA library construction.After library construction, the qualified libraries were sequenced using the Illumina NovaSeq 6000 using PE150 mode.Following an extracting and filtering quality control, we obtained high-quality, cleaned reads, and a follow-up analysis was then conducted (Table S1).All experiments were repeated three times with biological replicates.The statistical power of this experimental design, calculated in RNASeqPower is 0.96, based on a sequencing depth of 6 GB, CV of 0.4.We have uploaded the RNA-seq into the NCBI, the NCBI accession number is PRJNA971380.

DEGs and differential alternative splicing (DAS) analysis
We used FeatureCount (version 2.0.2) (Liao, Smyth & Shi, 2014) to quantify transcripts at the gene level.Differential expression analyses were performed with edgeR (version 3.3.3)according to the criteria of |log2 (FC)| > 1 and P value < 0.05.
Alternative splicing (AS) is the process by which different splice sites in precursor messenger RNA are selected to generate multiple mRNA isoforms, so AS is an important mechanism in creating proteome diversity and regulating gene expression in different tissues and developmental stages.To identify the number of different splicing events in TAO patients and controls, the software rMATS (version 4.0.2) was used (Shen et al., 2014), a new statistical method for robust and flexible detection of differential AS from replicate RNA-Seq data.Five main alternative splicing events, A3SS, A5SS, MXE, RI, and

Functional enrichment analysis
GO enrichment analyses for both the upregulated and downregulated genes were carried out using the R package topGO (The Gene Ontology Consortium, 2021) and the results were visualized using the REVIGO tool (http://revigo.irb.hr)(Supek et al., 2011).KEGG Orthology Based Annotation System (KOBAS) v3.0 (Bu et al., 2021) was used to perform the functional enrichment analysis.GSEA was carried out using the R package 'clusterProfiler' (Yu et al., 2012).The results are indicated in the appropriate figure legend and text.

The protein-protein interaction (PPI) network and hub gene identification
Construction of a PPI network was conducted using STRING (https://string-db.org/).We uploaded DEGs to STRING and obtained high-resolution bitmaps.By calculating the degree of connectivity, the hub genes in the PPI network were identified via cytoHubba, which is a plugin in Cytoscape software (version v3.9.1) (Shannon et al., 2003).

RNA quantification
Total RNA was extracted using the RNAiso Plus (Takara, Kyoto, Japan), according to the manufacturer's instructions.Final RNA pellets were resuspended in nuclease-free H2O and then determine the purity and concentration by measuring the optical density at 260 nm and 280 nm (NanoDrop 2000c; Thermo Fisher Scientific, Waltham, MA, USA).Reverse transcription of total isolated RNA was performed using the PrimeScript RT master mix kit (Takara, Kyoto, Japan).Gene expression was measured by qRT-PCR.The data were analysed using the 2 − CT method and normalized to the endogenous control GAPDH mRNA (for humans), and the amount of target gene mRNA expression in each sample was expressed relative to that of the control.Primer sequences for qRT-PCR were designed using Primer Express Software (Thermo Fisher Scientific, Waltham, MA, USA; Table S2).

Statistical analysis
The results are expressed as the mean ± SD.Significance was established between the two groups using Student's t test (paired t test).Age was compared using the t -test, and gender was compared using chi-squared tests.The data were analysed using GraphPad Prism 5 statistical software (Prism v5.0; GraphPad Software, La Jolla, CA, USA).A P value < 0.05 was considered statistically significant.

DEGs in orbital adipose/connective tissue samples of TAO patients
Deep sequencing identified 183 DEGs with the conditions of |log2(FC)| > 1 and P value < 0.05 between the orbital adipose/connective tissues of TAO patients and control individuals.Among these, 114 genes were upregulated, and 69 genes were downregulated.
The fragments per kilobase million (FPKM) value of mRNAs shows that there is no

DAS gene analysis
Alternative splicing (AS) refers to the process of selectively removing or retaining exons/introns during the initial transcription of DNA into RNA and further processing into mature mRNA, resulting in multiple transcripts of a gene.To learn the potential AS of TAO patients, five main types of AS events were analyzed using rMATS, including exon skipping (SE), intron retention (RI), alternative 5 splice site (A5SS), alternative 3 splice site (A3SS), and mutually exclusive exons (MXE) (Fig. 2A).We selected the DAS genes with a threshold of P value < 0.01.The numbers of A3SS, A5SS, MXE, RI, and SE events were 65, 57, 22, 18, and 477, respectively.SE was the most prevalent AS event in TAO patients, whereas RI was the least prevalent (Figs.2B, 2C).This data suggests that an abnormal splicing process leads to specific splicing isoforms, which may have a close relationship with the occurrence and development of TAO.

Enrichment analyses of DEGs
To explore the functions of DEGs, functional enrichment analysis was performed on DEGs by linking them with biological phenomena and their underlying mechanisms.

GO annotation analyses
GO analysis is a common useful method for large-scale functional enrichment research, which can significantly distribute DEGs into the biological process (BP), molecular function  2.
In the GO terms of TAO samples, inflammation response was the main BP category, including inflammatory response, regulation of inflammatory response, acute inflammatory response, regulation of acute inflammatory response, and myeloid leukocyte migration (Fig. 3A).This suggests that the pathogenesis of TAO is closely related to the aberrant activation of inflammatory responses, which play a key role in the activation of orbital adipogenesis.The MF category was abundant in glycosaminoglycan binding, G protein-coupled receptor binding, signaling receptor binding, and extracellular matrix structural constituent (Fig. 3B).In addition, CC mainly displayed extracellular region, extracellular space, and cell surface (Fig. 3C).

KEGG pathway enrichment analyses
The KEGG database is a widely used database to systematically analyze the metabolic pathways of gene products in cells and the functions of these gene products.It can help us study genes and expression information as a whole network.By analyzing the signaling pathway of DEGs, we can understand the significantly changed metabolic pathway in the state of TAO, which is important for exploring the pathogenesis of the disease.
KEGG analysis showed that 142 pathways were significantly enriched.The top 20 enriched pathways are shown in Fig. 3D.The represented pathways were ECM-receptor interaction, PI3K-Akt signaling pathway, cell adhesion molecules, cytokine-cytokine receptor interaction, and focal adhesion.

GSEA
GSEA is a promising, widely used software package that derives gene sets to determine different biological functions between two groups.By GSEA, we identified that cytokinecytokine receptor interaction, cytokine-cytokine receptor interaction, NF-kappa B signaling pathway, rheumatoid arthritis, TNF signaling pathway, and viral protein interaction with cytokine and cytokine receptor were the top five enriched pathways (Fig. 3E).In summary, the biological processes from the enriched GO terms, KEGG pathways, and GSEA for the DEGs were mainly involved in the regulation of inflammatory response, glycosaminoglycan binding and hyaluronic acid binding.

The protein-protein interaction (PPI) network and hub gene
To better understand the molecular mechanism of TAO, we visualized the importance of the relationship between proteins of DEGs using Cytoscape software (Table S3, Fig. 4B).Moreover, we identified the top 10 HUB genes according to node degree among these target genes via the Cytoscape plug-in cytoHubba (Fig. 4C).Collectively, these results suggest that the core proteins CXCL8, Toll-like receptor-2 (TLR2), CCL4 and angiotensinogen (AGT) in the PPI network may be involved in the regulation of TAO pathogenesis.

Histology and inflammation in the orbital adipose/connective tissues of TAO patients and control individuals
H&E staining showed the morphology of orbital adipose tissue, and consistently indicated an increased level of the inflammatory cell infiltration (black arrows) in TAO patients compared with the control individuals (Fig. 6A).Meanwhile, we immunohistochemically stained sections for the CD45, a protein expressed on all leukocytes, and found that CD45 expression (black arrows) also increased in the TAO patients compared with controls (Fig. 6B).To identify and quantitate macrophages within adipose tissue, we detected the expression of F4/80 antigen, a marker specific for mature macrophages in the orbital adipose tissues.Indeed, the TAO group had significantly increased amounts of F4/80-positive macrophages, compared with the control groups.
Besides, in the orbital muscle tissues, the inflammatory markers, CD45 and ICAM1 also increased surrounding the myofibrils in patients with TAO compared with the control groups (Figs.S3A and S3B).Moreover, there was a potent increase in the expression of fibrotic proteins, including α-SMA and FN (Figs.S3C and S3D), indicating the orbital fibrosis or myositis in individuals with TAO.
Collectively, our results showed that there were enhanced inflammatory responses in orbital adipose/connective tissue and increased levels of fibrosis in the extraocular muscles among TAO patients.

DISCUSSION
TAO is an autoimmune disease that affects orbital adipose tissue and extraocular muscles (Bartalena & Tanda, 2022).To date, the pathogenic mechanisms of TAO have not been clearly understood.Symptomatic treatments, such as hormone pulse therapy and orbital decompression, are currently limited for patients with TAO (Baeg et al., 2022).This issue emphasizes the importance of understanding the underlying mechanism(s) of and identifying therapeutic approaches for the prevention or treatment of TAO.In this study, we analysed the DEGs in orbital adipose/connective tissue from TAO patients and controls.
The symptoms of TAO are mainly caused by the inflammation in the orbital connective tissue, an increase in orbital volume due to enhanced adipogenesis and overproduction of glycosaminoglycans, and fibrosis of the extraocular muscles (Kahaly et al., 1992).
It has been reported that the inflammatory levels significantly upregulated in the adipose tissue and muscle of TAO patients (Carroll et al., 2013;Natesha et al., 1992;Khong et al., 2015).Huang et al. (2022) demonstrated that endoplasmic reticulum stress initiated by cholesterol metabolism may provoke adipose inflammation in TAO.Adipocyte-derived CP and AGT play a critical role in adipogenesis as well as inflammation (Carroll et al., 2013;Bednarek, Wysocki & Sowinski, 2004).Consistent with previous studies, we found elevated levels of CP and ATG in the adipose tissue and muscle of TAO patients.Existing data show that SERPINA3, an acute phase response protein, is involved in the pathogenesis of acute anterior uveitis, chronic obstructive pulmonary disease, Parkinson's disease, Alzheimer's disease, and coronary artery disease (Eidet et al., 2021;Li et al., 2021;Sánchez-Navarro et al., 2021).There is also literature supporting that SERPINA3 can be expressed to promote cell proliferation, migration, and expression of inflammatory cytokines by NF-κB signaling pathways (Liu et al., 2022).Consistently, SERPINA3 is also upregulated in both adipose tissue and muscle in TAO.Our research also combined RNA sequencing analysis with multiple validation experiments including qRT-PCR, H&E, immunohistochemistry and immunofluorescence analysis.H&E staining, CD45 and ICAM1 immunohistochemistry staining, and F4/80 immunofluorescence staining results showed the inflammatory responses potently increased in the orbital adipose/connective tissues of TAO patients, compared with the control groups (Figs. 6A and 6B,Figs. S3A and S3B).
In our study, RSPO1 was downregulated more significantly in orbital connective tissue than that in orbital fatty tissue.We speculated that this may be related to the fibrosis of the extraocular muscles.There is literature indicating that in other organs, such as the kidney, RSPO1 plays an important role in fibrogenesis, which may explain why the downward trend of RSPO1 is more pronounced in muscles (Su et al., 2021).
In our study, the results of GO molecular function analysis indicated that these DEGs were enriched in several terms, such as glycosaminoglycan binding, and extracellular matrix structural constituent.Wu et al. (2020), Wu et al. (2021a) and Wu et al. (2021b) indicated that several extracellular matrix related mRNAs (such as COL12A1, COL6A3) significantly reduced in TAO samples and closely related to the abnormal deposition of the extracellular matrix in orbital fat tissues in TAO patients (Liang et al., 2021).Additionally, GSEA and KEGG pathway enrichment analyses of the DEGs also showed marked enrichment of the NF-κB pathway, ECM-receptor interaction, cell adhesion molecules, and PI3K-Akt signaling pathway.During the pathogenesis of TAO, orbital fibroblasts are thought to interact with immunocompetent cells recruited to the orbit (Heufelder, 1995).They produce high amounts of glycosaminoglycans, particularly hyaluronan, which absorb water and lead to an increase in matrix volume (Smith et al., 1991).It has been documented previously that CD40-CD40 ligand interactions have important roles in the activation of human orbital fibroblasts (Cao et al., 1998).CD40L-provoked signaling pathways, including the NF-kappa B pathway and PI3K-Akt signaling pathway, result in the high expression of a variety of cytokines, such as VCAM-1, E-selective protein, IL-6, and other cytokines, in orbital fibroblasts of patients with TAO (Gillespie et al., 2012;Hwang et al., 2009).Fibroblasts are reported to be responsible for the secretion of collagen, release of extracellular matrix, and participation in inflammatory responses (Smith & Janssen, 2019).This functional characterization is further substantiated by α-SMA and FN immunofluorescent staining results (Figs.S3C and S3D).
As a previous study reported, there may exist alterations in the composition of the intestinal microbiota among patients, who suffered from severe and active TAO (Mori, Nakagawa & Ozaki, 2012).We found that pathway analyses highlighted the enrichment of highly expressed genes in the intestinal immune network for IgA production.In a separate investigation, Shi et al. (2019) found that two gut microbiotas (s_Prevotella_copri and f_Prevotellaceae) showed a significant correlation with TRAb.This suggests that intestinal symbiotic microorganisms may influence extraintestinal immune responses through the mucosal immune response induced by IgA antibodies, and they may render tolerance to self-antigens incompetent, such as TRAb, which can stimulate orbital and periorbital tissues and constitutes an independent risk factor for GO (Pianta et al., 2017;Seo & Sanchez Robledo 2018).
As with all transcriptomic analyses, there are limitations to this study.With the use of human tissue, there is heterogeneity in the patient's genetic background and other characteristics, such as age, gender, and CAS, which likely affect the disease.As such, we removed the influence of smoking on our results as much as possible, which has a strong and consistent association with TAO (Bartalena et al., 1989).One notable limitation lies in the relatively small sample size employed in our study, which consequently limits the statistical power.Additionally, while we selected genes that we believed were most important to the pathogenic mechanisms of TAO, it is imperative to acknowledge the presence of numerous other DEGs and pathways presented in these results that could be important and contribute to TAO.

CONCLUSIONS
Our transcriptome analysis identified 183 DEGs between TAOs and normal orbit tissues.Through an integrated bioinformatics analysis and verification of the DEGs, we identified several key candidate genes and enriched pathways that may aid the search for biomarkers and therapeutic targets of TAO.However, further molecular biology experiments are required to validate the findings of this study.

Figure 1
Figure 1 The differentially expressed genes were analyzed from RNA sequencing data.(A) Volcano plot of different genes in control or TAO orbital fat.FC, fold change; DEGs, differentially expressed genes.(B) Hierarchical clustering heatmap showing gene expression differences.Full-size DOI: 10.7717/peerj.16569/fig-1

Figure 2
Figure 2 Analysis of differential alternative splicing (AS) genes and distribution of the five main AS events.(A) Schematic diagrams of the mechanisms of the five main AS events.(B) Venn diagram of the detected genes undergoing the five AS events and overlap of these genes.SE, exon skipping; RI, intron retention; A5SS, alternative 5 splice site; A3SS, alternative 3 splice site; MXE, mutually exclusive exons.(C) Distribution of differential AS events based on a threshold of P < 0.01.Full-size DOI: 10.7717/peerj.16569/fig-2

Figure 3
Figure 3 The most significantly enriched GO terms and KEGG pathway analysis relevant to up-and downregulated genes.(A) BP term of GO enrichment analysis, *p < 0.001.BP, biological process.(B) MF term of GO enrichment analysis, *p < 0.05.MF: the molecular function.(C) CC term of GO enrichment analysis, *p < 0.05.CC: cellular component.(D) KEGG pathway analysis showing pathways that are enriched in the TAO group.(E) Gene cluster enrichment analysis (GSEA) revealed a significant enrichment of the first five pathways in TAO patients.Full-size DOI: 10.7717/peerj.16569/fig-3

Figure 4
Figure 4 The Venn diagram and the top hub genes identified in the protein-protein interaction (PPI) networks.(A) The Venn diagram shows the differentially expressed gene identification in the two gene expression profile datasets.(B) PPI network of differentially expressed genes.(C) Identification of the top 10 hub genes.Full-size DOI: 10.7717/peerj.16569/fig-4

Figure 5 Figure 6
Figure 5 Validation of the expression levels of mRNAs in the TAO groups and control groups.(A and B) The mRNA expression levels in orbital adipose tissue as verified by qRT-PCR.(C and D) Expression levels of mRNAs in orbital muscle tissues as verified by qRT-PCR.The results are presented as the means ± SDs; n = 4, * p < 0.05, and **p < 0.01 for each pair of groups indicated.Full-size DOI: 10.7717/peerj.16569/fig-5

Table 2 The top GO terms of DEGs between TAO and control samples.
The top 10 BP terms, MF terms and the most significantly CC terms of DEGs between TAO and control samples.